Dynamic systems with more than one dependent variable

Predator-Prey Models

Predator-Prey models

A simple approach that assumes prey grow exponentially, with a fixed intrinsic growth rate

Predator-Prey Model

Analogs

As with diffusion, the basic form/ideas in this model can be applied elsewhere

Differential equations for a Predator-Prey Model

\(\frac{\partial prey}{\partial t} = r_{prey} * prey - \alpha * prey * pred\)

\(\frac{\partial pred}{\partial t} = eff * \alpha * pred * prey - mort * pred\)

Predictor Prey - Implementation in R

Example implementation

source("../R/lotvmod.R")
lotvmod
## function (t, pop, pars) 
## {
##     with(as.list(c(pars, pop)), {
##         dprey = rprey * prey - alpha * prey * pred
##         dpred = eff * alpha * prey * pred - pmort * pred
##         return(list(c(dprey, dpred)))
##     })
## }
# note the use of with
# initial conditions
currpop = c(prey=10, pred=1)

# time points to see results
days = seq(from=1, to=100, by=1)

# set parameters
pars = c(rprey=0.5, alpha=0.3, eff=0.2,pmort=0.2, K=100)

# run the model
res = ode(func=lotvmod, y=currpop, times=days, parms=pars)

# graph the results
head(res)
##      time      prey     pred
## [1,]    1 10.000000 1.000000
## [2,]    2 11.320956 1.564997
## [3,]    3 10.244569 2.482924
## [4,]    4  6.914805 3.420960
## [5,]    5  3.777091 3.836373
## [6,]    6  1.987468 3.710744
# rearrange for easy plotting
ressimple = as.data.frame(res) %>% gather(key="animal", value="pop",-time)
p1=ggplot(ressimple, aes(time, pop, col=animal))+geom_line()

p1

Dynamics of two variables

Relationship between populations

p2=ggplot(as.data.frame(res), aes(pred, prey))+geom_point()
 p2

ggarrange(p1,p2)

# try other parameters to bring relative size of preditors higher

Other illustations

Preditor and Prey with Carrying Capacity

How would you code that?

source("../R/lotvmodK.R")
lotvmodK
## function (t, pop, pars) 
## {
##     with(as.list(c(pars, pop)), {
##         dprey = rprey * (1 - prey/K) * prey - alpha * prey * 
##             pred
##         dpred = eff * alpha * prey * pred - pmort * pred
##         return(list(c(dprey, dpred)))
##     })
## }
# initial conditions
currpop=c(prey=1, pred=1)

# set parameter list
pars = c(rprey=0.1, alpha=0.6, eff=0.8,pmort=0.4, K=20)

# times when you want to evaluate
days = seq(from=1,to=500)

# run our differential equation solver
res = ode(func=lotvmodK, y=currpop, times=days, parms=pars)

# estract the results
ressimple = as.data.frame(res) %>% gather(key="species", value="pop",-time)

# graph both populations over time
p1=ggplot(ressimple, aes(time, pop, col=species))+geom_line()
p1

# also look at relationships between preditor and prey population and use color for time 
# I will remove the legend here to make it easier to see 
p2 = ggplot(as.data.frame(res), aes(pred, prey, col=as.factor(round(time/10))))+geom_point()+theme(legend.position = "none")
p2

ggarrange(p1,p2)

# try with different parameter sets, can you create one where populations are stable - less cycling?

Sensitivity analysis

Consider pred-prey BUT what will be the output - if we want to ‘quantify sensitivity’ useful to look at a single value or set of value

For example Max Prey Pop Min Prey Pop

Sensitivity Analysis steps

source("../R/lotvmodK.R")
library(lhs)
library(epiR)
## Loading required package: survival
## Package epiR 2.0.46 is loaded
## Type help(epi.about) for summary information
## Type browseVignettes(package = 'epiR') to learn how to use epiR for applied epidemiological analyses
## 
# create parameter sets...
factors = c("rprey","K","alpha","eff","pmort")
nsets=500

# create a latin hyper cube
sens_pp = randomLHS(nsets, length(factors))
colnames(sens_pp)=factors

# refine using sampling distributions for each parameter
sens_pp[,"rprey"]= qunif(sens_pp[,"rprey"], min=0.01, max=0.3)
sens_pp[,"K"]= qunif(sens_pp[,"K"], min=10, max=200)
sens_pp[,"alpha"]= qunif(sens_pp[,"alpha"], min=0.1, max=0.4)
sens_pp[,"eff"]= qnorm(sens_pp[,"eff"], mean=0.3, sd=0.01)
sens_pp[,"pmort"]= qunif(sens_pp[,"pmort"], min=0.05, max=0.45)

# lets create a metric and wrapper function as we did for our population models

# first our metrics 
# lets say we  want the maximum and minimum  of both predictor and prey

compute_metrics = function(result) {
  maxprey = max(result$prey)
  maxpred = max(result$pred)
  minprey = min(result$prey)
  minpred = min(result$pred)
return(list(maxprey=maxprey, minprey=minprey, maxpred=maxpred, minpred=minpred))}

# build a wrapper function


p_wrapper = function(rprey,alpha, eff, pmort, K, currpop, days, func) {
    parms = list(rprey=rprey, alpha=alpha, eff=eff, pmort=pmort, K=K)
    result = ode(y=currpop, times=days, func=func, parms=parms) 
    colnames(result)=c("time","prey","pred")
  # get metrics
  metrics=compute_metrics(as.data.frame(result))
  return(metrics)
}


# run our model for all parameters and extract the results
currpop=c(prey=1, pred=1)
days = seq(from=1,to=500)
allresults = as.data.frame(sens_pp) %>% pmap(p_wrapper, currpop=currpop, days=days, func=lotvmodK)
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 379.35, R2 = 2.44102e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 379.35, R2 = 2.44102e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 379.35, R2 = 2.44102e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 379.35, R2 = 2.44102e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 379.35, R2 = 2.00841e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 379.35, R2 = 2.00841e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 379.35, R2 = 2.00841e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 379.35, R2 = 2.00841e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 379.35, R2 = 1.65247e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 379.35, R2 = 1.65247e-14
##  
## DLSODA-  Above warning has been issued I1 times.  
##      It will not be issued again for this problem.
## In above message, I1 = 10
##  
## DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
##       taken on this call before reaching TOUT     
## In above message, I1 = 5000
##  
## In above message, R1 = 379.35
## 
## Warning in lsoda(y, times, func, parms, ...): an excessive amount of work (>
## maxsteps ) was done, but integration was not successful - increase maxsteps
## Warning in lsoda(y, times, func, parms, ...): Returning early. Results are
## accurate, as far as they go
# take results back to unlisted form
allres = allresults %>% map_dfr(`[`,c("maxprey","minprey","maxpred","minpred"))


# range of response across parameter uncertainty
allresl = allres %>% gather(key="metric",value="pop")
ggplot(allresl, aes(metric, pop))+geom_boxplot()

# dealing with different scales
ggplot(allresl, aes(metric, pop, col=metric))+geom_boxplot()+facet_wrap(~metric, scales="free")

# plot cummulative densities

ggplot(allresl, aes(pop, col=metric))+stat_ecdf(geom="line")+facet_wrap(~metric, scales="free")

# compute PRCCs using epi.prcc
# lets do first for maxpred

epi.prcc(cbind.data.frame(sens_pp, allres$maxpred))
##            est       lower       upper test.statistic  df       p.value
## 1  0.535188659  0.46081666  0.60956065    14.13846752 498  2.131286e-38
## 2  0.446927930  0.36816815  0.52568771    11.14904975 498  6.358944e-26
## 3 -0.825948596 -0.87558190 -0.77631529   -32.69524914 498 4.732649e-126
## 4  0.004423757 -0.08361744  0.09246496     0.09872116 498  9.213994e-01
## 5  0.783353606  0.72862875  0.83807847    28.12406243 498 6.589430e-105
# try minprey
epi.prcc(cbind.data.frame(sens_pp, allres$minprey))
##           est       lower       upper test.statistic  df       p.value
## 1  0.82349307  0.77354450  0.87344165      32.392293 498 1.111315e-124
## 2 -0.05047556 -0.13840540  0.03745427      -1.127846 498  2.599280e-01
## 3 -0.78936918 -0.84341985 -0.73531850     -28.693524 498 1.363173e-107
## 4  0.02744035 -0.06056856  0.11544925       0.612587 498  5.404290e-01
## 5  0.65080575  0.58396033  0.71765118      19.128659 498  1.473946e-61

Stability

How do we think about stablity?

Populations don’t change when derivatives are zero!

What conditions lead to BOTH derivatives being zero

Stability For lotvmod

Make dprey and dpred equal to 0 and rearrange

Stability For lotvmodK

Make dprey and dpred equal to 0 and rearrange

Try setting you initial conditions close to these values and see what happens

# set parameter list
pars = data.frame(rprey=0.1, alpha=0.6, eff=0.8,pmort=0.4, K=200)

# now lets try initial conditions that will be stable
preyi = with(pars, pmort/(eff*alpha)) 
predi = with(pars, rprey/alpha*(1-preyi/K)) 
# times when you want to evaluate
days = seq(from=1,to=500)

# lets first see what happens when we start with 1 of each
currpop=c(prey=1, pred=1)
# run our differential equation solver
res = ode(func=lotvmodK, y=currpop, times=days, parms=pars)
# extract the results
res_smallstart = as.data.frame(res) %>% gather(key="animal", value="pop",-time)
# graph both populations over time
p1=ggplot(res_smallstart, aes(time, pop, col=animal))+geom_line()
p1

# lets first see what happens when we start our estimates of stable populations
stablepop = c(prey=preyi, pred=predi)
res = ode(func=lotvmodK, y=stablepop, times=days, parms=pars)
# estract the results
res_stablestart = as.data.frame(res) %>% gather(key="animal", value="pop",-time)
# graph both populations over time
p2=ggplot(res_stablestart, aes(time, pop, col=animal))+geom_line()
p2

# of course in this case these are not very realistic populations, so it is unlikely that you would get to something stable given parameters - but if they were different you might find a stable and realistic population

# try

And just for fun

Lets look at a Lorenz System with its interesting dynamics

# lorenze
source("../R/lorenz.R")
pars = list(a=10,b=28,c=8/3)
res = ode(func=lorenz, c(x=0.1,y=0,z=0), times=seq(0,50,by=0.01), parms=pars)

ggplot(as.data.frame(res), aes(x,y, col=time))+geom_point()

ggplot(as.data.frame(res), aes(x,z, col=time))+geom_point()

ggplot(as.data.frame(res), aes(y,z, col=time))+geom_point()

ressimple = as.data.frame(res) %>% gather(key="var", value="value",-time)
ggplot(ressimple, aes(time, value, col=var))+geom_line()

# try with different initial conditions
pars = list(a=15,b=28,c=8/4)
res = ode(func=lorenz, c(x=0.3,y=5,z=10), times=seq(0,50,by=0.01), parms=pars)

ggplot(as.data.frame(res), aes(x,y, col=time))+geom_point()+scale_colour_gradientn(colours = terrain.colors(10))

ggplot(as.data.frame(res), aes(x,z, col=time))+geom_point()+scale_colour_gradientn(colours = terrain.colors(10))

ggplot(as.data.frame(res), aes(y,z, col=time))+geom_point()+scale_colour_gradientn(colours = terrain.colors(10))

ressimple = as.data.frame(res) %>% gather(key="var", value="value",-time)
ggplot(ressimple, aes(time, value, col=var))+geom_line()